function [Delta_RCL, dDelta_RCL] = rcnl_delta( ...
    share_main,         ...
    Pi,            ...
    rho,           ...
    xnonlin_jt,          ...
    dfull_jt          ...
)

    J = size(xnonlin_jt,1);
    nmarket = size(xnonlin_jt,3);

    Delta_RCL         = zeros(J,1,nmarket);
    dDelta_RCL         = zeros(J,4,nmarket);

    for sel = 1:nmarket
        
        share_main_temp = share_main(:,:,sel);
        non_missing_index = share_main_temp ~= 0;
    
        outshr = (1-sum(share_main_temp(non_missing_index,:),1));
        logodds = log(share_main_temp(non_missing_index,:)) - log(outshr);
        expmval_mat = [tempname(), '.mat'];
        expmval0 = exp(logodds);
        save(expmval_mat, 'expmval0');
        
        cdindex = splitapply(@(s) s(end), (1:size(share_main_temp(non_missing_index,:), 1))', 1);
    
        [delta, ~] = rcnl_meanval( ...
            Pi, ...
            rho, ...
            share_main_temp(non_missing_index,:), ...
            xnonlin_jt(non_missing_index,:,sel), ...
            cdindex, ...
            dfull_jt(non_missing_index,:,sel), ...
            expmval_mat, ...
            1, ...
            1e-6 ...
        );

        % calculate ddelta 
        ddelta = rcnl_ddelta( ...
            Pi, ...
            rho, ... 
            delta, ...
            xnonlin_jt(non_missing_index,:,sel), ...
            cdindex, ...
            dfull_jt(non_missing_index,:,sel), ...
            1 ...
        );
        
        Delta_RCL(non_missing_index,:,sel) = delta;     
        dDelta_RCL(non_missing_index,:,sel) = ddelta;     
        
    end
end
